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ABSTRACT 



Aims. To determine whether or not gas-phase chemical models with homogeneous and time-independent physical conditions explain the many 
observed molecular abundances in astrophysical sources, it is crucial to estimate the uncertainties in the calculated abundances and compare 
them with the observed abundances and their uncertainties. Non linear amplification of the error and bifurcation may limit the applicability of 
i-C chemical models. Here we study such effects on dense cloud chemistry. 
Or Methods. Using a previously studied approach to uncertainties based on the representation of rate coefficient errors as log normal distributions, 
q . we attempted to apply our approach using as input a variety of different elemental abundances from those studied previously. In this approach, 
!— I all rate coefficients are varied randomly within their log normal (Gaussian) distribution, and the time-dependent chemistry calculated anew 
C/j , many times so as to obtain good statistics for the uncertainties in the calculated abundances. 

■ Results. Starting with so-called "high-metal" elemental abundances, we found bimodal rather than Gaussian like distributions for the 
abundances of many species and traced these strange distributions to an extreme sensitivity of the system to changes in the ratio of the cosmic 
• i— I , ray ionization rate £ H e for He and that for molecular hydrogen £ H2 • The sensitivity can be so extreme as to cause a region of bistability, 
1 which was subsequently found to be more extensive for another choice of elemental abundances. To the best of our knowledge, the bistable 
5-H | solutions found in this way are the same as found previously by other authors, but it is best to think of the ratio fae/ify as a control parameter 
. . . perpendicular to the "standard" control parameter f/« H - 



Key words. Astrochemistry - ISM: abundances - ISM: clouds - ISM: molecules 



1. Introduction 

The equations describing the chemical evolution of quiescent 
cores (also known as molecular clouds) are highly non-linear 
and can result in extreme sensitivity to the initial conditions 
or to some of the parameters. For example, the elemental 
ratio of carbon to oxygen is well known to be an important 
param eter for dense cloud chemistry jTerzieva & H erbst 
Il998h . Another well-known consequence of this non-linearity 
is the presence of bistable regions for some ranges of model 
parameters (temperature, density, comic -ray ionization rate, 
rate coefficients and elemental abundances). The phenomenon 
of bistability is characterized by two stable solutions for 
chemical abundances at steady-state for the same set of 
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parameters. Bistability in dark cloud co nditions, i.e. at 
low temperature, was first discovered by iLe Bourlot et alJ 



authors jLe Bourlot et al 
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1995 


Lee et al. 


1998 




3 ineau Des Forets & Roueff 


2000; 


ICharnlev & Markwickl 120031 iBoeer & Sternberd 12006). The 



two solutions of this bistability are characterized by a large 
difference in the C/CO and O/O2 ratios and in the ionization 
fraction. For this last reason, they were named the H igh 
and L ow Ionization Phases (HIP and LIP). ILe Bourlot et alJ 
dl995l) showed that the region of bistability can be mapped 
out by variations in density, temperature, cosmic -ray ioniza- 
tion rate, and elemental d e pletion s within specific ranges. 
Pineau Des Forets & Roueff (2000) later showed that even 
variations in rate coefficients can le ad to the two solutions. 
Recently, iBoger & Sternberg! J2OO6) identified the chemical 
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Table 1. High- and low- metal elemental abundances with re- 
spect to H2. 



Element 


High-Metal 


Low-Metal 






(Intermediate) 


He 


2.8 x lfr 1 


2.8 x 10- 1 


N 


4.28 x 10~ 5 


4.28 x 10~ 5 


O 


3.52 x 1(T 4 


3.52 x 10~ 4 


C + 


1.46 x 10~ 4 


1.46 x 10~ 4 


s + 


1.6 x 1(T 5 


1.6 x 10~ 7 






(2.0 x lO" 6 ) 


Si + 


1.6 x 10- 6 


1.6 x 10~ 8 


Fe + 


6.0 x 1(T 7 


6.0 x 10-" 


Na + 


4.0 x 10~ 7 


4.0 x 1Q- 9 


Mg + 


1.4 x 10- 7 


1.4 x 10- 8 


P + 


6.0 x 10~ 7 


6.0 x 10-' 


cr 


8.0 x 10~ 7 


8.0 x 1Q- 9 



Note: The intermediate-metal case differs from the low-metal case 
by the S + abundance only, as is indicated in brackets in the table. 



mechanisms of this phenomenon as a cycle involving the H3 , 
O2, and S + species. 

In two previous papers, we presented a Monte Carlo 
method to compute the theoretical error of abundances in 
chemical models d ue to unc ertainties in rate c oefficients and 
physical conditions (Wakelam et al. 2005, 2006). The errors in 
the rate coefficients are defined by a log-normal distribution. As 
a consequence, the resulting abundances at steady-state usually 
follow a Gaussian distribution ( Wak elam et alJ 120051) . When 
applying this method, we found that the chemical abundances 
are very sensitive to the ratio between the ionization rates of 
helium fee and molecular hydrogen £h 2 caused by cosmic -ray 
particles. Characterized by a variation of several orders of mag- 
nitude in some abundances for small variations of ^He/fife, this 
sensitivity can lead to bimodal rather than Gaussian distribu- 
tions for the abundances. For certain ranges of the parameters, 
it can even lead to bistability, with the ratio <THe/<TH 2 appearing 
to be a control parameter. 

This paper, which describes the sensitivity and some of its 
consequences, is organized as follows. In section 2, we briefly 
present the chemical model and the uncertainty method used 
for this analysis. In a somewhat unorthodox way, we have de- 
cided to present, in section 3, the manner in which the sensitiv- 
ity was found while studying the use of so-called "high-metal" 
elemental abundances. Section 4 shows the influence of the el- 
emental abundances and the chemical network used on the sen- 
sitivity. In section 5, we present and discuss the phenomenon of 
bistability as obtained with the variation of fHe/fH 2 - In section 
6, the chemical reactions involved in this sensitivity/bistability 
are elucidated. Finally, the last section presents our conclu- 
sions. 

2. Chemical model and method of uncertainty 

We used a gas-phase time-depen dent chemical mod el with the 
osu.2003 1 network reported by Smi th et alJ J2004 (4233 re- 

1 http://www.physics.ohio-state.edu/~eric/research_files/ 



actions, 421 species and 12 elements). The model computes 
the evolution of the species for a fixed temperature and den- 
sity. The initial conditions are all atomic except for molecu- 
lar hydrogen, w hile the chemical mod el and database are the 
same as used in Wakela m et alJ J2006I) . For this work, we will 
present the results using three different sets of elemental abun- 
dances: the high-, low- and intermediate-metal cases. The high- 
and low-metal elem ental abundances have been defined by 
iGraedel etalJ fl982) and are listed in Tabled The intermediate 
case has the same abundances as the low-metal one except that 
the amount of the element sulfur is raised to 2 x 10~ 6 compared 
with H2 (also given in Tabled- The other parameters are the 
typical ones for quiescent cores: a kinetic temperature of 10 K, 
an H2 density of 10 4 ctrT 3 , a visual extinction of 10 so that the 
photochemistry driven by the external UV photons does not oc- 
cur, and a fixed cosmic -ray ionization rate ( of 1 .3 x 10~ 17 s _1 . 

The method used in this analysis to derive uncertaintie s 
in abundance is described in detail in Wakela m et alJ (|Jo05). 
The uncertainties in rate coefficients derive from the UMIST 
database when available, but for most of the reactions, includ- 
ing the ionization rates of H2 and He, we have assumed a factor 
of 2 uncertainty. Note that the uncertainty in ionization rates 
is not an uncertainty in the parameter but an uncertainty in 
the multiplicative factor that distinguishes the individual rates. 
Although ionization rates for H2 and He have been used in a 
large number of models, the curren t rates used in terms of £ 
derive back to the Ph. D. thesis of IB lack! Jl975i) . where the 
factors are 0.93 for the production of H+ (£°,) and 0.5 for the 
production of He + (£° e ), leading to a ratio of 0.54, if 

one neglects t he minor channels fo r H?. A n older ratio of unity 
was used bv iHerbst & Klemperen .19731) . New estimates, in- 
volving both direct ionization and secondary ionization by en- 
ergetic electrons, are sorely needed (Blac k, priva te communi- 
cation) based on the treatment of Dalgarn o et alJ (|l_999) for a 
mixture of atomic and molecular hydrogen and helium. 

The uncertainty method consists of generating /V new sets 
of rate coefficients by replacing each coefficient k, by a ran- 
dom value consistent with its uncertainty factor F,-. We as- 
sume a normal distribution of log k, with a standard deviation 
cr, = log Fj. We run the model for each set j, which produces, 
for each species, N values of the fractional abundances Xj(t) at 
a time t. For this work, we ignore the uncertainty in physical 
conditions and consider the uncertainty in the rate coefficients 
only. A total /V of 2000 different runs was made for each set 
of parameters studied; this number is large enough to achieve 
statistical significance (see Wakel am et all2 006). 

3. The sensitivity to £h 6 /£h 2 in the high-metal case 

The computation of the theoretical error due to the uncertainties 
in rate coefficients was first applied to the case of high-metal el- 
emental abundances. Figure^shows the results of the random 
variation of all the rate coefficients for the O2 abundance as a 
function of time. As can be seen from this figure, the distribu- 
tion of the abundance is spread into two peaks after 10 6 yr, as 
steady-state approaches, with a significant number of curves in 
each peak and some stable solutions between the two peaks. 
The bimodality of the distribution at a specific time (10 7 yr) is 
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Fig. 1. Density of probability of the O2 abundance as a function 
of time (left panel). The right panel shows the histogram of 
the abundance at 10 7 yr. The elemental abundances are for the 
high-metal case. 
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Fig. 2. Distributions (nomalized to the total number of runs) of 
the abundances of C, H2O, e~ and H3 at 10 7 yr (high-metal 
case). The solid lines refer to the variation of all the rate coef- 
ficients whereas the dotted lines represent the variation of ^ e 
and £n 2 only. The vertical dashed line represents the abundance 
computed with the standard values of the rate coefficients. 



also shown as a histogram in Figure Such bimodal distribu- 
tions are obtained for a large number of species (C, H2O, SO, 
CS etc) with the largest separation between the two peaks for 
atomic carbon (three orders of magnitude). Fig. |2 shows his- 
tograms for the abundance distributions of the species C, H2O, 
e~, and HJ at 10 7 yr. It can be seen that the H^ distribution 
does not show a bimodal profile but is Gaussian. The same is 
true for H2S and OH. Since the error in the rate coefficients fol- 




log(k[cm-S]) 



Fig. 3. Histograms of the rate coefficients of the following reac- 
tions (upper panel): He + CRP — > He + + e~ (£ne, left box), H2 
+ CRP -> H+ + e" (£h 2 , middle box) and H+ + e -> H 2 + H 
(right box). The thick and thin lines represent the histograms of 
all the runs (hist 1) and of the runs giving the higher abundance 
of O2 (hist 2, see text) respectively. The lower panels show the 
ratio between hist 2 and hist 1 . 



lows a log-normal distribution, the bimodal shapes are due to 
non-linear effects. 

We obtained the bimodal distributions by varying the rate 
coefficients, we then decided to identify the reactions responsi- 
ble for them if at all possible. We proceeded as follows. Among 
the 2000 runs, we listed the sets of reactions that give one the 
higher of the two most probable abundances for O2. Then for 
each rate coefficient, we plotted the histograms for all runs and 
the histograms for those runs giving only the higher abundance 
of O2. In the upper panels of Fig. [5] pairs of histograms can be 
seen for the rate coefficient of the cosmic ray ionization of He, 
that of H2, and the dissociative recombination of H^ to form H2 
and H. The thick lines (histl) represent histograms for all runs 
while the thin lines (hist2) represent those ending up with the 
higher O2 abundance. For reactions not responsible for the bi- 
modal distributions, there should be no difference between the 
pairs of histograms except that the total number of runs should 
be lower for the second one. For most of the reactions, we in- 



deed found no difference, as for the reaction Hi 



H 2 + 



H in Fig. |3] In fact, there were only two exceptions: the direct 
ionization rate coefficients of He and H2 by cosmic rays (fne 
and £h 2 ). and their marked effect can be seen in Fig. [3] The 
high abundance of O2 tends to be formed with lower values of 
^He and higher values of £n 2 , as can be seen most clearly in the 
lower panels of Fig. [3] in which the ratios of the two histograms 
are plotted. It is then convenient to consider the sensitivity of 
the abundance of O2 and other species to the ratio £ke/£H 2 rather 
than to their absolute values. To confirm their importance, we 
randomly varied only ^He and £h 2 , keeping the other rate coeffi- 
cients constant, and we obtained the same bimodal distribution 
as with the fully random method at steady state. The dispersion 
of the abundances in this case is then almost as large as when 
we varied all the rate coefficients, as can be seen in the dotted 
lines in Fig. [5] 

The sensitivity of the chemical abundances to the ratio 
^He/^H 2 can be explored in another way. In Fig. @] (left panel), 
we show the steady-state abundance of O2, computed by ran- 
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Fig. 4. Abundance of Cb as a function of ^He/^H 2 a t 10 7 yr. The 
upper and middle panels represent the abundance computed 
with the random distribution of all the rate coefficients and of 
^He and only, respectively. The solid lines in the lower plot 
have been computed using one value of ^He for each line and 
varying (u 2 > while the dotted lines have been computed us- 
ing one value of £h 2 f° r each nne and varying ^He- The vertical 
dashed line indicates the standard value of ■ 
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Fig. 5. Abundances of C, e~, H2O and H3 as a function of 
^He/^H 2 and for five different values of £ne (solid lines). The 
vertical dashed line indicates the standard value of £[L/£^ ■ 



the inflection point, confiming that the choice of the ionization 
ratio as a parameter of the sensitivity is a reasonable one. 

The abundance of O2 decreases by two orders of magni- 
tude (from 1(T 5 to 1(T 7 ) when ( H Jiu 2 goes from 0.3 to 0.8, 
and the inflection point of the curve occurs at the standard 
value of Qjg^ used in osu.2003 (0.5/0.93). In Fig.0we show 
some other examples of this sensitivity by using five differ- 
ent values of gn e , as in the solid lines of the right panel of 
Fig. |4] The abundance of atomic carbon has an almost equal 
sensitivity albeit opposite to that of O2, while the H2O abun- 
dance shows a similar profile to O2 except for a local maxi- 
mum at ^He/^H, around 0.2. Although the electronic abundance 
increases towards higher £He/^H 2 , it varies only a factor of 2.5 
between ^ H e/^H 2 =0.1 and ^He/^H 2 =3. The abundance of the H3 
ion, which does not show any bimodal distribution, is linear 
with ^He/^k 2 - The dispersion of the abundances, for a given ra- 
tio ^He/fH 2 . due to the variation of ^He by a factor of 2 around 
the standard value £jL, depends on the species: it is smaller for 
H2O, O2 and e~ than for C and H3 , for instance, but remains 
anyway below 0.5 dex. 



domly varying all the rate coefficients as a function of ^He/^H, ■ 
For comparison, the abundance computed with the random 
variation of ^He and £h 2 only is shown in the middle plot of the 
figure. The smaller vertical dispersion is caused by not varying 
the other rate coefficients. Finally and for more clarity, we de- 
pict, in the lower panel, the calculated O2 abundance as a func- 
tion of ^He/^H, obtained by calculations in which the ionization 
rates are not varied randomly. Rather, the solid lines show val- 
ues of the steady-state O2 abundance computed for 5 values of 
£ He (between 2£^ e and ^ /2) with £ H , varied to account for the 
range of the abscissa between 0.05 and 5. In the same panel, 
dotted lines show the O2 abundances computed for 5 values of 
with ^He varied to account for the range of the abscissa. The 
abundances computed by each method are similar especially at 



4. Sensitivity for different model parameters 

4.1. Other choices of elemental abundances 

The calculated sensitivity depends on the elemental abun- 
dances used for the modeling. In Fig. [6] we show the C, H2O, 
e~ and H3 fractional abundances as a function of £ke/fH 2 for 
the low-, intermediate-, and high-metal elemental abundances 
(see Sect.|2). The sensitivity plots, as well as most others later 
in the paper, are obtained by fixing the helium ionization rate 
at its standard value so that no dispersion is seen. The H2 den- 
sity remains at 10 4 cirT 3 and the temperature at 10 K. Although 
the sensitivity exists for the three cases, the inflection point is 
shifted to higher values of £He/£k 2 f° r the low-metal case. For 
the intermediate case, on the other hand, the sensitivity seems 
to be stronger, leading to an abrupt jump in the abundance of 
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Fig. 6. Steady-state abundances of C, e~, H2O and as 
functions of £ke/£H 2 f° r three different elemental abundances. 
The triangles, diamonds and crosses represent the low-, 
intermediate-, and high-metal cases. The vertical dashed line 
indicates the standard value of ^L/fg • The H2 density is 
10 4 ctrT 3 and the temperature 10 K. The value of ^He is fixed at 
its standard value while £h 2 is varied. 
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Fig. 7. Distribution (nomalized with respect to the total number 
of runs) of the steady-state abundance of atomic carbon for the 
three sets of elemental abundances. The vertical dashed line 
represents the abundance computed with the standard values of 
the rate coefficients. 



C for ^He/^H 2 -1-1- A similar abrupt change can be seen for the 
electronic abundance. This jump is actually a manifestation of 
bistability (see section|5|i. Note that for ^He/^H, < 1.1, the abun- 
dances for the intermediate-metal case tend to be close to the 
low-metal ones whereas for £nJ£n 2 > 1-1 they are similar to 
the high-metal ones. 

When one runs the uncertainty method with a factor of two 
uncertainty in £ne an d & 2 > the values of (bb/Chu are spread 
approximatively between 0.13 and 2.2. The distribution of 
the abundances computed with this method thus depends on 
the sensitivity within this range of values. The distribution 
of the steady-state abundances of C computed by variation 
of all the rate coefficients are shown in Fig. for the three 
sets of abundances. Because the inflection point for the low- 
metal case is shifted to £ H e/fH 2 =2, the uncertainty method gives 
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Fig. 8. Atomic carbon abundance at steady state as a function of 
£hc/&2 f° r the three elemental abundances (low-, intermediate- 
, and high-metal). In each plot, the symbols refer to different 
C/O ratios. The vertical dashed line indicates the standard value 



Gaussian-like distributions although slightly asymmetrical to- 
wards higher C abundances. This explains why we did not 
see the sensitivity effect in our previo us work on rate coeffi - 
cient uncertainties in dark clouds ("see IWakelam et aD E006). 
For high-metal elemental abundances, the point of inflection 
occurs close to the standard value of £ee/^H 2 =0.54 and the 
distribution for the C abundance is bimodal. When using the 
intermediate-metal abundances, the uncertainty method gives 
two distinct peaks with no solutions between them. The peak 
at higher abundance is the smaller since the point of inflection 
occurs at values of ^He/^H 2 =1-1 somewhat larger than the stan- 
dard. The lack of any model runs with steady-state abundances 
in between the two distributions stems from the sharpness of 
the curve around the inflection point and is a consequence of 
the bistability. 

In addition to its dependence on metallicity, the sensitivity 
to depends strongly on the C/O elemental abundance 

ratio. In Fig. [8] we show the steady-state abundance of atomic 
carbon as a function of fjiJfa computed for the three sets of 
elemental abundances, each with three different values of C/O 
- 0.2, 0.4 (the standard one), and 0.6 - obtained by maintain- 
ing a constant carbon abundance while varying the oxygen. The 
lower the C/O ratio is, the more the inflection point is shifted 
towards higher values of ^He/^Hi ■ As an example, the inflection 
point for the high-metal case occurs around 0.2 for C/O=0.6 
whereas it rises to 1 .5 for C/O=0.2. In both cases, the inflection 
point is too far from the standard value of £nd£n 2 to obtain a 
bimodal distribution. For the intermediate-metal case, the C/O 
ratio of 0.6 results in shifting the inflection point to the stan- 
dard value of £hJ(h 2 > so that the analogous panel to the one in 
Fig. would show two equal areas for the two peaks. The dis- 
persion of the carbon abundance caused by the variation of C/O 
is larger for values of (hJ(h 2 smaller than the inflection point 
while the abundances eventually form a single plateau after the 
jump. 

As discussed in Wakelam & Herbst (in preparation), the 
helium abundance used in the classic low- and high-metal el- 
emental abundances (see Section is quite high (0.28 com- 
pared with H2) compared with the more modern value of 
0.18 measured in t he interstellar medium (see, for example, 
iBaldwin et ai1ll99ll) . Naturally, the sensitivity will change if 
a lower value of He is used. Fig.[9]shows the distribution of the 
steady-state C abundance computed with a random variation 
of ^He and £h 2 and the high-metal case with the lower helium 



6 



V. Wakelam et al.: Extreme sensitivity of the dense cloud chemistry 




+ high metal 
O inter, metal 
A low metal 



j+ . O A 

\aaaa,a a 



-9 -8 -7 -6 -5 
log[Xl 



0.1 



1.0 

we* 



10.0 



Fig. 9. The left panel shows the distribution of the steady-state 
abundance of atomic carbon computed with a random variation 
of ^He and £h 2 f° r the high-metal elemental abundances but with 
a lower helium elemental abundance of 0.18 compared with 
H2. The right panel represents the C steady-state abundance as 
a function of ^He/^H, for the three elemental abundances with 
the low He abundance. The vertical dashed line indicates the C 
abundances computed with the standard and ^ on the left 
plot and the standard value of on the right plot. 




Fig. 11. Illustrative scheme of hysteresis surface demonstrat- 
ing bistability with two "orthogonal" control parameters, here 
ne, and ^ne/(n 2 )- The position of unstable solutions is indicated 
by dashed lines . Orthogonal cuts crossing the folding region 
produce hysteresis curves like the ones shown in Fig. 1121 and 
Fig. El 
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Fig. 10. Abundance of C as a function of £nd£n 2 computed with 
the high- (black lines) and low- (grey lines) metal abundances. 
The left and right panels show the results using the Srates (see 
text) and rate99 networks. The vertical dashed line indicates 
the standard value ^ e /^ 2 - 



abundance. In this case, we do not obtain a clear bimodal distri- 
bution but merely an asymmetrical profile, because the inflec- 
tion point is shifted towards higher values of fne/£H 2 > as shown 
in Fig. [5] and the random method barely reaches the sensitive 
portion of the curve. 

4.2. Chemical database 

To test the influence of the chemical network utilized, we con- 
structed the sensitivity curves at steady-state with two other 
databases - Srates and rate99- using both low-metal and high- 
metal abu ndances. The rate99 network is the UMIST list of 
reactions dLe Teuff et alJ EoOoT) prior to its current updating 
(see http://www.udfa.neti whereas Srates is the list of reactions 
used by Wakelam et al. J2004I) to study sulphur chemistry in 
hot cores. This latter network cont ains only five element s (H, 
He, C, O and S) but, as discussed in Wakelam et alJ (E004). the 
steady-state abundances computed with it are very close to the 
ones computed with the larger osu.2003 network at low tem- 



perature. The ionization rates for H2 and He are the same in all 
three databases. 

Fig- Ellis a sensitivity plot of the C abundance for five dif- 
ferent values of ^He between 2(n e and ^He/2. The sensitivity 
obtained with Srates is similar to that for osu.2003, which indi- 
cates that the absence of the other elements (Si, Fe, Na etc.) 
does not change the sensitivity drastically. On the contrary, 
the results with rate99 show a much different relation between 
^He/^Hi and the atomic carbon abundance, especially for the 
high-metal case, where the rate99 curves show (a) a signifi- 
cantly lower value for the point of inflection of the sensitivity, 
(b) a wide dispersion due to the value of £ne at low values of 
Cnd&h, an d ( c ) a muc h smaller jum p in the abundance. We 
have shown in Wakel am etailfcOOfil) that rate99 and osu.2003 
often have differences in rate coefficients larger than the esti- 
mated random errors for unstudied reactions of a factor of 2, 
whereas most of the rate coefficients in Srates are taken from 
the osu database. Since the three networks have the same ^He 
and £h 2 > the shape of the sensitivity can depend only on the 
values of other rate coefficients. 



5. Bistability 

In section 4, we saw that the variation of (hJ(h 2 leads to a sharp 
change in the abundances of some species for the intermediate- 
metal case (see Fig. [6). This jump is indeed a m anifestation of 
the same bistability discovered by Le Bou rlot et al.l dl993h . as 
can be determined by comparing the region of the bistability 
in parameter space. The parameter usually used to study this 
bistability is the ratio between the cosmic -ray ionization rate £ 
and the density «h but other model parameter s such as the tem- 
perature and rate coefficients have been u sed dLe Bourlot et alJ 
ll995UPineau Pes Forets & RoueflfcOOol . What we have found 
here is that ^ke/^H, is also a control parameter which explore 
another dimension of the bistability. If we imagine £ife/& 2 to 
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Fig. 12. Hysteresis curves of the steady-state abundance of 
atomic carbon as a function of (hJ(h 2 f° r three different ele- 
mental abundances (low-, high-, and intermediate-metal), and 
for five H2 densities between 10 2 and 10 6 cirT 3 . The vertical 
dashed line indicates the standard value of £jL/£^ ■ Grey zones 
indicate regions of bistability. 
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Fig. 13. Hysteresis curves of the electron steady-state abun- 
dance as a function of (left panel) and fke/fe. (right panel). 
Intermediate-metal elemental abundances are used for both 
plots. Five H2 densities between 10 2 and 10 6 cirT 3 are used for 
the right panel. The vertical dashed line indicates the standard 
value of ■ Grey zones indicate regions of bistability. 



be an orthogonal axis to £/«h, the bistability o flLe Bourlot et all 
Jl993h is the same as ours when £He/fH 2 is at its standard value. 
In Fig. ^2 we show an illustrative scheme of this bistability as 
controlled by the two parameters nn, and ^He/^H, ■ To make the 
bistability appear, we constructed hysteresis curves as a func- 
tion of ^He/^Ho for different H2 densities. These curves, shown 
in Fig.^lfor atomic carbon, were produced in two steps. Fist, 
we ran the model to compute the steady-state abundance us- 
ing one value of ^He/^H 2 > then we used this solution as the ini- 
tial condition and increased tjijtji 2 to compute the next point, 
shown as the crosses in Fig. [Q] and [D] In addition, the same 
principle but with decreasing £nJ%H 2 was utilized, with results 
shown as diamonds.. In order to vary £n e /£n 2 , we fixed the he- 
lium ionization rate and changed that of hydrogen. When two 



solutions exist for a range of £hJ£h^ values, a region of bista- 
bility exists. 

To explore the domains of bistability using the new control 
parameter, we constructed hysteresis curves using the three dif- 
ferent elemental abundances and five densities between 10 2 and 
10 6 crrT 3 . All the results are shown in Fig.^lfor atomic car- 
bon. Bistability exists for all the elemental abundances but the 
domain of ^ke/^Hj where the two solutions coexist are small. To 
show the bistability more clearly, a magnification of this region 
is done for n(H2)=10 4 cirT 3 and the intermediate-metal case. 
In both the low- and high-metal cases, bistability exists only 
at very low density (10 2 crrT 3 ) whereas in the intermediate- 
metal case, the bistability exists between 10 2 and 10 4 cirT 3 . 
The increase in density shifts and smoothes the inflection point, 
making the bistability eventually disappear. Note that this is 
the first ti me that th e bistabil ity has been fo und with the high- 
metal case JLee et alll998tlBoger & Sternbergfcooiil) . because 
it does not exist for the standard value of fHe/£k 2 ■ We have con- 
firmed that, for the high-metal case, using the standard value 
of (hJ(b 2 leads to no bistability with nH, as low as 10 cirT 3 . 
Using ^He/^H 2 = 0.16 however, bistability exists in a very small 
density range of 70-90 cirT 3 . 

A more focused look at the two control parameters is shown 
in Fig. where we plot the hysteresis curves of the elec- 
tron abundance as a function of each parameter for the case 
of intermediate-metal elemental abundances. The value of £ 
is set at the standard value. For the variation of £n e /£n 2 (right 
panel), we used five H2 densities between 10 2 and 10 6 cm -3 , 
as in the previous figure, while for the hysteresis curve as 
a function of nn 2 (left panel) the standard value of fne/fn 
(0.54) was used and the bistability found to exist in the density 
range 60-150 cm -3 . With the other control parameter and for 
nn 2 = 10 2 crrT 3 , the bistability exists for ^W^ks between 0.5 
and 0.7 and the domains of bistability using both parameters 
overlap. For higher densities, the bistable region using £ke/£H 2 
shifts towards higher values and so does not exist anymore in 
the hysteresis curve versus ne, . These hysteresis curves can be 
combined into a plot of the type shown in Fig. II II 

6. Chemical reactions involved in the sensitivity 

The sensitivity of the abundances to the ratio £hJ(h 2 is re- 
lated to the non-linearity of the system. A detailed mathemat- 
ical analysis of this phenomenon is in progress (Massacrier et 
al., in preparation). It is possible, however, to understand the 
difference, between the low and high £He/£n 2 phases qualita- 
tively from a chemical point of view. For this purpose, we have 
considered the main reactions of production and destruction of 
O, C, CO, S and SO, which are the main reservoirs of car- 
bon, oxygen and sulphur for the intermediate- and low-metal 
elemental compositions in which sulfur is relatively abundant, 
and/or are the species most influenced by ^He and £h 2 ■ Let us 
first consider low values of (ue/(u^ ■ The ionization of H2 leads 
to the formati on of OH through a well- known series of reac- 
tions (see. e.g. iBoger & Sternbergl2006l) . The radical OH then 
reacts with atomic oxygen to form O2. The most efficient way 
of destroying O2 is to form SO by reaction with S. Since the el- 
emental abundance of atomic sulphur is much less than that of 
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oxygen, large amounts of O2 remain in the gas phase. Atomic 
carbon reacts with O2 to form CO where it is stored since CO 
is not efficiently destroyed. For high ^ke/fifei on the other hand, 
the destruction of CO by He + produces C + and O more effi- 
ciently. Carbon ion then exchanges its charge with S so large 
abundances of C remain in the gas phase to deplete molecular 
oxygen. 

Since the bistability we see is an extreme manifestation of 
the sensitivity to ^He/^H 2 , the chemical processes involved in 
bistability may be related to those discussed above, although 
our contr ol parameter is dif ferent fro m those previously stud- 
ied.lPineau des Forets et a H dl992l) and lLe Bourlot et alJ Jl993l 
1995) identified the abundance ratio H + /HJ as being impor- 
tant for bistability, with one ion identified with the HIP and 
the other with the LIP solution. The abundance ratio H + /HJ 
depends on the abundance of electrons and as a consequence 
on the efficiency of the dissociative recombination of HJ with 
electrons, the rate of which was quite uncertain at the time. In 
our analysis however, the H3 ion shows little if any sensitivity 
to the control parameter, although the electron abundance does 
show a jump between the low and high ^He/^H, phases, which 
can then be thought of as an alogous in some sense to th e LIP 
and HIP solutions. Recently, Boger & Sternber made 
another detailed analysis of the chemical differences between 
the so-called HIP and LIP solutions, arguing that bistability is 
due to a loop including H^, O2 and S + . Although this expla- 
nation is related in some aspects to that given for the sensitiv- 
ity reported here, one must remember that our control parame- 
ter raises some new points, given the relation of the ionization 
rate of He to the formation of the important ion C + , which is a 
precursor for many species, and can charge exchange to form 
abund ant C. Moreover, S + , which is prominent in the explana- 
tion of iBoger & Sternberg ! J2006). is not the dominant charge 
carrier in our calculations. Note that the osu.2003 database in- 
cludes the recombination of the main ions with negative grains 
and this does not cancel the bistability, although such a cancel- 
lation was suggested by Boger & Sternberg ( 2006). It may well 
be that the H + /Ht abundance ratio, the H^ -02-S + cycle, and the 
^He/& 2 ratio comprise pieces of a larger and more complex jig- 
saw puzzle. Since bistability also appears to be a function of 
the reaction network utilized, secondary reactions may play a 
more important role than heretofore realized. This possibility 
has to be clarified through the building of a simplified reaction 
network that still shows the main characteristics of bistability. 

7. Observational relevance 

The sensitivity to ^ke/^H, depends highly on the model param- 
eters and, for this reason, it is not obvious if such effect can be 
detected in the interstellar medium. In an ideal situation where 
the cloud characteristics (at least the elemental abundances and 
the age) are well known and correspond to a distinct chemical 
regime, one could use the observations to constrain the ratio 
^He/^H^ assuming however that no local variation of this ratio is 
expected. In reality, the cloud characteristics are in general far 
from being well constrained, and an analysis of the uncertain- 
ties in the ionization rates of H2 and He is strongly required to 
understand the relevance of this sensitivity for the interstellar 



medium. Note that £nj£n 2 is n °t sensitive to the total intensity 
of the cosmic -ray flux. 

If the intermediate-metal elemental abundances reflect the 
cloud characteristics and if the ionization rate of He and H2 by 
cosmic -rays can vary by at least a factor of 2 within a cloud, 
the prediction of bistability over a wide range of densities sug- 
gests that strong chemical heterogeneities in dense clouds can 
be caused by this effect. In the absence of bistability an ex- 
treme sensitivity is predicted at most densities for low-metal 
and high-metal abundances and may cause some observed het- 
erogeneities in abundance. However, we have to consider the 
age of the cloud. The bifurcation occurs just after 10 5 yr, start- 
ing from our standard initial conditions, and is a maximum at 
steady-state (> 10 7 yr). As a consequence, the chemical hetero- 
geneity can be caused by this sensitivity only in evolved clouds 
O 10 5 - 10 6 yr). 

The large gas-phase abundances of H2O and O2 molecules 
predicted by dense cloud chemical models is a long standing 
problem. Using high values of fHe/& 2 » we obtain much smaller 
abunda nces for these species, closer to the upper limits in dense 
clouds (Bergi n & Snell2002t|Pagani et all2003l) . Large uncer- 
tainties in the £nJ(n 2 ra tio could then solve the problem re- 
lated to the modeling of water and molecular oxygen. One 
has, however, to consider the overall picture, and the better 
agreement with H2O and O2 should not worsen the agree- 
ment for the other species. A comparison between observations 
and m odeling s hould be done in a systematic way, as done by 
Wakelam et al. (2006). This was not the purpose of the present 
study, but we can already note, as an example, that the observed 
abundances in L134N (North peak) of O H and C 4 H (7.5 x 10~ 8 
and 10~ 9 respectively, Ohishi et al. 1992) are reproduced by the 
model when using the standard value of ^He/^H, but are under- 
estimated when using the higher ratios ^ke/fe that reproduce 
H2O and O2 abundances. 

Finally, the very low H2 density (100 cm -3 ), at which the 
bistability as a function of (uJ(n 2 1S restricted to with the low- 
and high-metal elemental abundances, is not relevant for clouds 
with the assumed visual extinction (A,,=10) since it would re- 
quire a size of 26 pc, which is much larger than typical cloud 
sizes dSnellll98lh . 

8. Conclusions 

Using a Monte Carlo uncertainty method to compute the theo- 
retical error of molecular abundances in chemical models, we 
found a strong sensitivity of the steady-state abundances to the 
ratio between the ionization rate coefficient of helium ^He and 
that of molecular hydrogen £h 2 by cosmic rays. This sensitivity 
to the ratio has the consequence of changing the abundance of 
some species by several orders of magnitude for small varia- 
tions of ^He/^Hi ■ Lower values of ^He/^H 2 are characterized by 
large abundances of smaller molecules such as O2, H2O, SO, 
etc. whereas high values of ^He/^H, lead to large abundances 
of atoms. If the sensitivity is strong enough, it can result in 
the phenomenon of bistability, which we located at low densi- 
ties for two sets of elemental abundances - the low- and high- 
metal cases, and over a wide range of densities for the case 
of the so-called intermediate-metal abundances, in which sul- 
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fur plays a particularly prominent role. The bistability has been 
determined to be the same as studied earlier by a variety of au- 
thors with other control parameters; it is perhaps simplest to 
consider our control parameter - <fHe/£H 2 - as orthogonal to one 
of the others - the ratio of the overall ionization rate parameter 
£ to the density n - so that the bistability can be imagined to 
exist in a region in three-dimensional space rather than in the 
customary two-dimensional hysteresis picture. 

It is interesting to speculate that the sensitivity detected 
here via our uncertainty analysis is a necessary but not suffi- 
cient condition for bistability because only the most extreme 
sensitivity, in which the bimodal distributions detected have no 
solutions in between them, leads to bistability. Unfortunately, 
it is not guaranteed that our type of uncertainty analysis can 
always locate a bistable region since finding an extreme sensi- 
tivity may require large uncertainties in some parameters, such 
as the gas density, to take us from our standard parametric val- 
ues to the bistable region. 

Several questions arise when considering the dramatic ef- 
fects of small ^He/^H, variations on the gas-phase chemistry. 
First, even if one assumes a given flux and energy distribution 
of the cosmic rays, the determination of <5W<Th^ suffers some 
significant uncertainties, with the role of secondary electrons 
particularly important. Besides these difficulties, one might 
also have to consider variations of the ratio £He/<fH 2 m space, 
such as between one cloud and another, or between the cen- 
ter and the edges of a same cloud, and in time. Changes in the 
ratio imply important variati ons in the spec tral distribution of 
the cosmic rays (Glass gold & Langerlll973l) . Such variations 
might be produced by the attenuation of the cosmic rays inside 
a cloud or the proximity of supernovae in star forming clouds. 
Whether such variations are expected or not and how they could 
affect the chemistry should be considered in the future. Finally, 
a knowledge of the helium elemental abundance appears to be 
crucial since the sensitivity region depends on this parameter. 
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